*
* $Id$
*
* $Log: ghesig.F,v $
* Revision 1.1.1.1  2002/06/16 15:18:40  hristov
* Separate distribution  of Geant3
*
* Revision 1.1.1.1  1999/05/18 15:55:19  fca
* AliRoot sources
*
* Revision 1.1.1.1  1995/10/24 10:21:15  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 29/03/94  15.41.37  by  S.Giani
*-- Author :
      FUNCTION GHESIG(PPART,EKIN,AVER,A,Z,W,KK,DENS,QCOR,LPART)
C
C *** CALCULATION OF THE PROBABILITIES FOR (IN)ELASTIC INTERACTIONS ***
C *** OF STABLE PARTICLES ON PROTON AND NEUTRON                     ***
C *** NVE 07-APR-1988 ***
C
C CALLED BY : GPGHEI
C ORIGIN : F.CARMINATI, H.FESEFELDT (ROUTINE INTACT 06-OCT-1987)
C
C *** IPART DENOTES THE GHEISHA PARTICLE INDEX ***
C
C CONVENTION :
C
C   PARTICLE                 IPART
C   ------------------------------
C   GAMMA                    1
C   NEUTRINO                 2
C   POSITRON                 3
C   ELECTRON                 4
C   MUON +                   5
C   MUON -                   6
C   PION +                   7
C   PION 0                   8
C   PION -                   9
C   KAON +                  10
C   KAON 0 S                11
C   KAON 0 L                12
C   KAON -                  13
C   PROTON                  14
C   PROTON BAR              15
C   NEUTRON                 16
C   NEUTRON BAR             17
C   LAMBDA                  18
C   LAMBDA BAR              19
C   SIGMA +                 20
C   SIGMA 0                 21
C   SIGMA -                 22
C   SIGMA + BAR             23
C   SIGMA 0 BAR             24
C   SIGMA - BAR             25
C   XSI 0                   26
C   XSI -                   27
C   XSI 0 BAR               28
C   XSI - BAR               29
C   DEUTERON                30
C   TRITON                  31
C   ALPHA                   32
C   OMEGA -                 33
C   OMEGA - BAR             34
C   NEW PARTICLES           35
C
C   HE3                     36 (49)
C   ANTIDEUTERON            37 (76)
C   ANTITRITON              38 (77)
C   ANTIHE3                 39 (78)
C   ANTIALPHA               40 (79)
C
#include "geant321/gcunit.inc"
#include "geant321/gconsp.inc"
C
C --- FOR CROSS-SECTION INFORMATION SEE "PCSDATA" ---
C
#include "geant321/gsecti.inc"
C --- GHEISHA COMMONS ---
#include "geant321/s_result.inc"
#include "geant321/s_prntfl.inc"
C
      DIMENSION A(KK),Z(KK),W(KK)
C
      DIMENSION ALPHA(40),ALPHAC(41),IPART2(7),CSA(4)
      DIMENSION PARTEL(40),PARTIN(40),ICORR(40),INTRC(40)
C
      PARAMETER (ONETHR=1./3.)
      SAVE ALPHA,ALPHAC,PARTEL,PARTIN,CSA,IPART2,ICORR,INTRC
C
#include "geant321/s_csdim.inc"
#include "geant321/pcodim.inc"
#include "geant321/s_csdat.inc"
#include "geant321/pcodat.inc"
C
      DATA ALPHA    / 6*0.7,
     +                0.75 ,0.75 ,0.75 ,
     +                0.76,0.76 ,0.76 ,0.76 ,
     +                0.685,0.63 ,0.685,0.63,0.685,0.63,
     +                3*0.685,3*0.63,2*0.685,2*0.63,
     +                3*0.7,0.685,0.63,0.7,5*0.7/
      DATA ALPHAC    /1.2,1.2,1.2,1.15,0.90,0.91,0.98,1.06,1.10,1.11,
     +                1.10,1.08,1.05,1.01,0.985,0.962,0.945,0.932,
     +                0.925,0.920,0.920,0.921,0.922,0.923,0.928,0.931,
     +                0.940,0.945,0.950,0.955,0.958,0.962,0.965,0.976,
     +                0.982,0.988,0.992,1.010,1.020,1.030,1.040/
      DATA PARTEL/6*0.,29*1.,5*1./
      DATA PARTIN/6*0.,1.00,0.00,1.05,1.20,1.35,1.30,1.20,1.00,1.30,
     +            1.00,1.30,1.00,1.30,1.00,1.00,1.00,1.30,1.30,1.30,
     +            1.00,1.00,1.30,1.30,1.00,1.,1.,1.,1.3,1.,5*1./
      DATA ICORR /14*1, 0, 1, 0, 1, 0, 3*1, 3*0, 2*1, 2*0, 4*1, 2*0,
     +            5*1/
      DATA INTRC /6*0, 1, 0, 12*1, 0, 2*2, 0, 10*1, 0, 5*1/
C
C CROSS SECTIONS ON NUCLEUS ARE KNOWN ONLY FOR PIONS AND PROTONS.
C THE GENERAL LAW SIGMA(A)=1.25*SIGMA(TOT,PROTON)*A**ALPHA IS VALID
C ONLY FOR MOMENTA > 2 GEV.THE PARAMETRIZATION DONE HERE GIVES ONLY
C A BEHAVIOUR AVERAGED OVER MOMENTA AND PARTICLE TYPES.
C FOR A DETECTOR WITH ONLY A FEW MATERIALS IT'S OF COURSE MUCHBETTER
C TO USE TABLES OF THE MEASURED CROSS SECTIONS .
C FOR ELEMENTS WITH THE FOLLOWING ATOMIC NUMBERS MEASURED CROSS
C SECTIONS ARE AVAILABLE (SEE "PCSDATA").
C
C                 H   AL     CU     PB
      DATA  CSA  /1. ,27.00 ,63.54 ,207.19 /
      DATA IPART2/9,8,7,11,10,13,12/
C
C
C --- Initialise GHESIG and switch to GHEISHA particle code ---
      GHESIG=0.0
      IF(LPART.GT.79)GO TO 160
      IPART=KIPART(LPART)
C
C --- No interaction for gammas, neutrinos, electrons, positrons, muons,
C --- neutral pions, neutral sigmas and antisigmas and new particles.
      IF(INTRC(IPART).EQ. 0) GO TO 160
      P=PPART
      EK=EKIN
C
C --- Initialise the cross-sections with 0.0 ---
      DO 10  K=1,KK
         AIEL(K)=0.0
         AIIN(K)=0.0
         AIFI(K)=0.0
         AICA(K)=0.0
   10 CONTINUE
C
      IF(((IPART.GE.30).AND.(IPART.LE.32)).OR.
     $   ((IPART.GE.36).AND.(IPART.LE.40))) THEN
C
C --- light (anti)nuclei reaction cross section parametrization ---
C
         DO 20 K=1,KK
            AIIN(K) = SIGNUC(LPART,EK,A(K))
   20    CONTINUE
         IF (NPRT(9)) THEN
            WRITE(CHMAIL,10000)
            CALL GMAIL(0,0)
         END IF
C
      ELSE IF ((IPART .EQ. 16) .AND. (EK .LE. 0.0327)) THEN
C
C --- Use tables for low energy neutrons ---
C --- get energy bin ---
         JE2=17
         DO 30 J=2,17
            IF (EK .LT. ELAB(J)) THEN
               JE2=J
               GO TO 40
            END IF
   30    CONTINUE
C
   40    JE1=JE2-1
         EKX=MAX(EK,1.0E-9)
         DELAB=ELAB(JE2)-ELAB(JE1)
         DO 70 K=1,KK
C
C --- Get A bin ---
            JA2=15
            DO 50 J=2,15
               IF (A(K) .LT. CNLWAT(J)) THEN
                  JA2=J
                  GO TO 60
               END IF
   50       CONTINUE
C
   60       JA1=JA2-1
            DNLWAT=CNLWAT(JA2)-CNLWAT(JA1)
C
C --- Use linear interpolation or extrapolation by Y=RCE*X+RCA*X+B ---
C
C --- Elastic cross section ---
C --- E interpolation or extrapolation at JA1 ---
            DY=CNLWEL(JA1,JE2)-CNLWEL(JA1,JE1)
            RCE=DY/DELAB
C --- A interpolation or extrapolation at JE1 ---
            DY=CNLWEL(JA2,JE1)-CNLWEL(JA1,JE1)
            RCA=DY/DNLWAT
            B=CNLWEL(JA1,JE1)-RCE*ELAB(JE1)-RCA*CNLWAT(JA1)
            AIEL(K)=RCE*EK+RCA*A(K)+B
C
C --- Inelastic cross section ---
C --- E interpolation or extrapolation at JA1 ---
            DY=CNLWIN(JA1,JE2)-CNLWIN(JA1,JE1)
            RCE=DY/DELAB
C --- A interpolation or extrapolation at JE1 ---
            DY=CNLWIN(JA2,JE1)-CNLWIN(JA1,JE1)
            RCA=DY/DNLWAT
            B=CNLWIN(JA1,JE1)-RCE*ELAB(JE1)-RCA*CNLWAT(JA1)
            AIIN(K)=RCE*EK+RCA*A(K)+B
C
            IZNO=Z(K)+0.01
            AICA(K)=11.12*CSCAP(IZNO)/(EKX*1.0E6)**0.577
   70    CONTINUE
         IF (NPRT(9)) THEN
            WRITE(CHMAIL,10100)
            CALL GMAIL(0,0)
         END IF
      ELSE
C
C --- Use parametrization of cross section data for all other cases ---
C
         IF (NPRT(9)) THEN
            WRITE(CHMAIL,10200)
            CALL GMAIL(0,0)
         END IF
C
C --- Get momentum bin ---
         J=40
         DO 80 I=2,41
            IF (P .LT. PLAB(I)) THEN
               J=I-1
               GO TO 90
            END IF
   80    CONTINUE
C
C --- Start with  cross sections for scattering on free protons ---
C --- use linear interpolation or extrapolation by Y=RC*X+B     ---
   90    DX=PLAB(J+1)-PLAB(J)
C --- Elastic cross section ---
         DY=CSEL(IPART,J+1)-CSEL(IPART,J)
         RC=DY/DX
         B=CSEL(IPART,J)-RC*PLAB(J)
         AIELIN=RC*P+B
C --- Inelastic cross section ---
         DY=CSIN(IPART,J+1)-CSIN(IPART,J)
         RC=DY/DX
         B=CSIN(IPART,J)-RC*PLAB(J)
         AIININ=RC*P+B
         ALPH=ALPHA(IPART)
         IF(IPART.LT.14) THEN
            DY=ALPHAC(J+1)-ALPHAC(J)
            RC=DY/DX
            B=ALPHAC(J)-RC*PLAB(J)
            CORFAC=RC*P+B
            ALPH=ALPH*CORFAC
C
            IPART3=IPART2(IPART-6)
C
C --- Elastic cross section ---
            DY=CSEL(IPART3,J+1)-CSEL(IPART3,J)
            RC=DY/DX
            B=CSEL(IPART3,J)-RC*PLAB(J)
            XSECEL=RC*P+B
C --- Inelastic cross section ---
            DY=CSIN(IPART3,J+1)-CSIN(IPART3,J)
            RC=DY/DX
            B=CSIN(IPART3,J)-RC*PLAB(J)
            XSECIN=RC*P+B
C
         END IF
C
         DO 100 K=1,KK
            AIEL(K)=AIELIN
            AIIN(K)=AIININ
C
            IF (A(K) .GE. 1.5) THEN
C
C --- A-dependence from parametrization ---
               CREL=1.0
               CRIN=1.0
C --- Get medium bin  1=Hydr.  2=Al  3=Cu  4=Pb ---
               I=3
               IF (A(K) .LT. 50.0) I=2
               IF (A(K) .GT. 100.0) I=4
               IF ((IPART .EQ. 14) .OR. (IPART .EQ. 16)) THEN
C
C --- Protons and neutrons ---
C
C --- Elastic cross section ---
                  DY=CSPNEL(I-1,J+1)-CSPNEL(I-1,J)
                  RC=DY/DX
                  B=CSPNEL(I-1,J)-RC*PLAB(J)
                  XSECEL=RC*P+B
C --- Inelastic cross section ---
                  DY=CSPNIN(I-1,J+1)-CSPNIN(I-1,J)
                  RC=DY/DX
                  B=CSPNIN(I-1,J)-RC*PLAB(J)
                  XSECIN=RC*P+B
                  IF (AIEL(K) .GE. 0.001) CREL=XSECEL/(0.36*AIEL(K)*
     +            CSA(I)**1.17)
                  AITOT=AIEL(K)+AIIN(K)
                  IF (AITOT .GE. 0.001) CRIN=XSECIN/(AITOT*CSA(I)**
     +            ALPH)
C
               ELSEIF (IPART .LT. 15) THEN
C
C --- Calculate correction factors from values on Al,Cu,Pb for all ---
C --- mesons use linear interpolation or extrapolation by Y=RC*X+B ---
C --- Note that data is only available for pions and protons
                  WGCH=0.5
                  IF (A(K) .LT. 20.0) WGCH=0.5+0.5*EXP(-(A(K)-1.0))
                  AIEL(K)=WGCH*AIEL(K)+(1.0-WGCH)*XSECEL
                  AIIN(K)=WGCH*AIIN(K)+(1.0-WGCH)*XSECIN
C
C --- This section not for kaons ---
                  IF (IPART .LT. 10) THEN
C
C --- Elastic cross section ---
                     DY=CSPIEL(I-1,J+1)-CSPIEL(I-1,J)
                     RC=DY/DX
                     B=CSPIEL(I-1,J)-RC*PLAB(J)
                     XSPIEL=RC*P+B
C --- Inelastic cross section ---
                     DY=CSPIIN(I-1,J+1)-CSPIIN(I-1,J)
                     RC=DY/DX
                     B=CSPIIN(I-1,J)-RC*PLAB(J)
                     XSPIIN=RC*P+B
C
                     IF (AIEL(K) .GE. 0.001) CREL=XSPIEL/(0.36* AIEL(K)
     +               *CSA(I)**1.17)
                     AITOT=AIEL(K)+AIIN(K)
                     IF (AITOT .GE. 0.001) CRIN=XSPIIN/(AITOT*CSA(I)
     +               **ALPH)
                  END IF
               END IF
               AIIN(K)=CRIN*(AIIN(K)+AIEL(K))*A(K)**ALPH
               AIEL(K)=CREL*0.36*AIEL(K)*A(K)**1.17
               AIEL(K)=AIEL(K)*PARTEL(IPART)
               AIIN(K)=AIIN(K)*PARTIN(IPART)
            END IF
  100    CONTINUE

C  ---------------------------------------------
C   Fix for low energy antiprotons from
C   Astroparticle Physics Volume 6,
C   Issues 3-4, April 1997, Pages 379-386
         IF(IPART.EQ.15) THEN
            DO 105 K=1,KK
               AIIN(K)= SIGAPA(EK,A(K))
  105       CONTINUE
         END IF
C  ----------------------------------------

C
C --- Fission cross sections ---
C --- A-dependence given by  sigma(3 MeV)=-67.0+38.7*Z**(4/3)/A ---
      END IF
      I=21
      DO 110 II=1,21
         IF (EK .LT. EKFISS(II)) THEN
            I=II
            GO TO 120
         END IF
  110 CONTINUE
C
  120 CONTINUE
      DO 130 K=1,KK
C
C --- No fission for materials with A < 230 ---
         IF (A(K) .GE. 230.) THEN
C
C --- Only data for U(233), U(235) and Pu(239) for Ek .le. 0.01 GeV ---
            J=4
            IF (EK .LE. 0.01) THEN
C
C --- Distinguish U(233), U(235), Pu(239) and U(238) and rest by J=1,4
               IF ((Z(K) .EQ. 92.0).AND.(ABS(A(K)-233.0).LT.0.5))THEN
                  J=1
               ELSEIF((Z(K).EQ.92.0).AND.(ABS(A(K)-235.0).LT.0.5))THEN
                  J=2
               ELSEIF((Z(K).EQ.94.0).AND.(ABS(A(K)-239.0).LT.0.5))THEN
                  J=3
               END IF
            END IF
            IF(J.EQ.4) THEN
C
               Z43BA=Z(K)**(4.0*ONETHR)/A(K)
               Z43BA=MAX(-67.0+38.7*Z43BA,0.)
            ELSE
               Z43BA=1.
            ENDIF
C
C --- Energy dependence taken from U(238) ---
C --- approximated as step-function ---
            AIFI(K)=CSFISS(J,I)*Z43BA
         END IF
  130 CONTINUE
C
C --- Corrections for compounds ---
C --- These corrections should only be applied to inorganic scintill. ---
C --- apply the correction only if user selected it within GEANT ---
      IF (QCOR .GT. 0.0.AND.KK.GT.1) THEN
C --- Do not apply corrections for anti-baryons ---
         IF (ICORR(IPART).EQ.1) THEN
C
C --- ACC40 between 0.3 and 0.5 for Pi+ and Pi0, 0.2 for other mesons ---
            ACC40=0.325
C           IF (IPART .GE. 9) ACC40=0.20
C --- ACC40 = 0.15 for baryons ---
            IF (IPART .GE. 14) ACC40=0.15
C --- ACCA=0.08 for all pions, 0.02 for all other particles ---
            ACCA= 0.08
            IF (IPART .GE. 10) ACCA=0.02
C
            ACCB=0.32*(ACC40-ACCA)
            ACC=ACCA-ACCB*LOG(EK)
            IF (ACC .GT. 0.5) ACC=0.5
            IF (ACC .GT. 0.0) THEN
C
               CAVER=AVER**ACC
               CAVER=1.+(CAVER-1.)*QCOR
               DO 140 K=1,KK
                  AIEL(K)=AIEL(K)*CAVER
                  AIIN(K)=AIIN(K)*CAVER
                  AIFI(K)=AIFI(K)*CAVER
                  AICA(K)=AICA(K)*CAVER
  140          CONTINUE
            END IF
         END IF
      END IF
C
C --- Calculate interaction probability ---
C
C --- Correction factor for high (P > 100 GeV/C) energies ---
      CORH=1.0
C --- Assume a LOG(P) dependence of the correction factor with values ---
C P = 100 GeV/C  ==> CORH = 1.
C P =   1 TeV/C  ==> CORH = 1.25
C     DX=LOG(1000.)-LOG(100.)
C     DY=1.25-1.
C     RC=DY/DX
C     B=1.-RC*LOG(100.)
C     CORH=RC*LOG(P)+B
      IF (P .GT. 100.) CORH=0.1085736156*LOG(P)+0.5
      ALAM=0.0
      DO 150 K=1,KK
         AFACT=AVO*1E-3*DENS*W(K)/A(K)
         AIEL(K)=MAX(CORH*AIEL(K)*AFACT,0.)
         AIIN(K)=MAX(CORH*AIIN(K)*AFACT,0.)
         AIFI(K)=MAX(CORH*AIFI(K)*AFACT,0.)
         AICA(K)=MAX(CORH*AICA(K)*AFACT,0.)
C
         ALAM=ALAM+AIEL(K)+AIIN(K)+AIFI(K)+AICA(K)
  150 CONTINUE
C
C --- Pass the interaction probability to GEANT ---
      GHESIG=ALAM
C
      GO TO 999
C
C --- Printout of skipped particles in case of interface debug ---
  160 IF (NPRT(9)) THEN
         WRITE(CHMAIL,10300) IPART
         CALL GMAIL(0,0)
      END IF
10000 FORMAT(' *GHESIG* GEOM X-SECT. FOR INEL. SCAT. OF D,T AND ALPHA')
10100 FORMAT(' *GHESIG* X-SECT. FROM LOW ENERGY NEUTRON TABLES')
10200 FORMAT(' *GHESIG* X-SECT. FROM PARAMETRIZATION OF DATA')
10300 FORMAT('0*GHESIG* GHEISHA PARTICLE ',I3,' SKIPPED')
  999 END
